Nucleation and Bulk Crystallization in Binary Phase Field Theory 
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We present a phase field theory for binary crystal nucleation. In the one- component limit, quan- 
titative agreement is achieved with computer simulations (Lennard- Jones system) and experiments 
(ice-water system) using model parameters evaluated from the free energy and thickness of the in- 
terface. The critical undercoolings predicted for Cu-Ni alloys accord with the measurements, and 
indicate homogeneous nucleation. The Kolmogorov exponents deduced for dendritic solidification 
and for "soft-impingement" of particles via diffusion fields are consistent with experiment. 
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Understanding alloy solidification is of vast practical 
and theoretical importance. While the directional geom- 
etry in which the solidification front propagates from a 
cool surface towards the interior of a hot melt is under- 
stood fairly well, less is known of equiaxial solidification 
that takes place in the interior of the melt. The latter 
plays a central role in processes such as alloy casting, hi- 
bernation of biological tissues, hail formation, and crys- 
tallization of proteins and glasses. The least understood 
stage of these processes is nucleation, during which seeds 
of the crystalline phase appear via thermal fluctuations. 
Since the physical interface thickness is comparable to 
the typical size of critical fluctuations that are able to 
grow to macroscopic sizes, these fluctuations are nearly 
all interface. Accordingly, the diffuse interface models 
lead to a considerably more accurate description of nu- 
cleation than those based on a sharp interface [|. 

The phase field theory, a recent diffuse interface ap- 
proach, emerged as a powerful tool for describing com- 
plex solidification patterns such as dendritic, eutectic, 
and peritectic growth morphologies ||. It is of interest 
to extend this model to nucleation and post-nucleation 
growth including diffusion controlled "soft-impingement" 
of growing crystalline particles, expected to be responsi- 
ble for the unusual transformation kinetics recently seen 
during the formation of nanocrystalline materials Q. 

In this Letter, we develop a phase field theory for crys- 
tal nucleation and growth, and apply it to current prob- 
lems of unary and binary equiaxial solidification. 

Our starting point is the free energy functional 
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developed along the lines described in ]|, |(|. Here </> 
and c are the phase and concentration fields, /(0, c) — 
WTg(4>) + [1 - P (</))] fs + P(4>)fL is the local free energy 
density, W = (1 — c)Wa + cWb the free energy scale, 
the quartic function g(<p) = (j) 2 (l — qb) 2 /4 that emerges 
from density functional theory |?J ensures the double- 
well form of /, while the function P{(j>) = (/> 3 (10 — lbqb + 



6</> 2 ) switches on and off the solid and liquid contributions 
fs,L, taken from the ideal solution model. (A and B refer 
to the constituents.) 

For binary alloys the model contains three parameters 
e, Wa and Wb that reduce to two (e and W) in the 
one-component limit. They can be fixed if the respective 
interface free energy 7, melting point X/, and interface 
thickness S are known j8|. Such information is available 
for the Lennard- Jones, ice- water, and Cu-Ni systems [[)), 
offering a quantitative test of our approach. 

Relying on the isothermal approximation the time 
evolution is described by Langevin equations dqb/dt = 
-M^SF/S^) + C/> and dc/dt = V[M C V(SF/Sc)} + Cvc 
where (SF/ Sx) stands for the functional derivatives (x = 
(/>, c), M x are mobilities, while and (vc are appropriate 
noises added to RHS to mimic thermal fluctuations. The 
dimensionless form of these equations is obtained 
by measuring length and time in units £ and £ 2 /Dl, 
t = t£ 2 /D L , r = r/f. Here £ and D L are the charac- 
teristic length scale and the diffusion coefficient in the 
liquid, while the quantities with tilde are dimensionless. 

The critical fluctuation (nucleus) is a non-trivial time- 
independent solution of the governing equations. For 
spherical symmetry (a reasonable assumption), the phase 
field equation reduces to V 2 (/> = A/i(0, c)/(e 2 T). Here 
A/i(0, c) = WTg'(ab) + [(1 - c ) A/ A + cA/ B ]P'(0) is the 
local chemical potential difference relative to the initial 
liquid, prime stands for differentiation with respect to 
the argument, the local concentration is related to the 
phase-field as c(0) = c oc e _2y /(l — + c^e -27 ), where 
y = v(W B - W A )g(0)/R + v(Af B - A/ A )[P(0) - 1]/BT, 
while Afi are the volumetric free energy differences be- 
tween the pure liquid and solid phases. Solving these 
equations numerically under boundary conditions that 
prescribe bulk liquid properties far from the fluctuations 
(0 — > 1, and c — > Coo for r — » 00), and zero field-gradients 
at the center, one obtains the free energy of critical fluctu- 
ation as W* = F — Fo. Here F is obtained by numerically 
evaluating Eq. (1) after having the time-independent so- 
lutions inserted, while Fq is the free energy of the initial 
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FIG. 1: Nucleation in the modified Lennard- Jones system: 
(a) Radial phase field profiles for critical fluctuations at sev- 
eral temperatures, (b) Relative interfacial free energy vs. re- 
duced temperature, (c) Comparison of nucleation rates pre- 
dicted by the phase field theory (PFT), the classical sharp 
interface theory (CNT) and computer simulations (squares) 
[|l5| . Short-dashed lines show the limits of nucleation rate al- 
lowed by the error of the interfacial free energy, e = 137.07/c 
and a = 3.383 A are taken for Ar. Below 61.7 K, the simula- 
tions increasingly underestimate the true nucleation rate due 
to an unknown equilibration period caused by quenching the 
liquid to the nucleation temperature Jli|]. 



liquid. This is compared with W* = (16n/3)tf/ Af 2 
from the sharp interface "droplet" model of the classical 
nucleation theory p| , where 7/ = j(Tf). 

The homogeneous nucleation rate is calculated as J = 
Joexp {— W*/kT}, where the nucleation prefactor Jo of 
the classical kinetic approach is used which proved 
consistent with experiments |TT] . 

To study the soft-impingement problem we introduce 
a non-conservative orientation field 0, which is random in 
the liquid, and has a constant value between and 1 in 
the crystal that determines crystal orientation in the lab- 
oratory frame. By this, we capture the feature that the 
short-range order in the solid and liquid are usually sim- 
ilar, with the obvious difference that the building units 
have a uniform orientation in the crystal, while their ori- 
entation fluctuates in the liquid. Following [|l3|, we as- 
sume that the grain boundary energy acts in the solid 
and is proportional to |V0|. We realize this by adding 
fori — .M|V0| to fs, where coefficient M. is assumed to 
be independent of c. The respective equation of motion 
has the form 80/ dt = -M 9 (SF/S0) + to, yielding 
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where to = te$P{4>), and M e = M e ,s + P{4>){M e , L - 
Mo y s), while subscripts S and L indicate the values for 
the bulk liquid and solid phases. When (j) < 1, Eq. 
(2) switches in orientational ordering, and chooses the 
value of that survives as the orientation of the particle, 



which then serves as the direction relative to which the 
anisotropy of 7 (= 70 {1 + socos[m(i!) — 6)]} for m-fold 
symmetry) is measured. A similar model has been suc- 
cessfully applied for describing grain boundary dynamics 
g| 0- Unlike previous work, in our approach the ori- 
entation field is coupled to the phase field, and extends 
to the liquid phase where crystallographic orientation de- 
velops from orientational fluctuations. While our model 
incorporates grain boundary dynamics, our primary in- 
terest is solidification, and Mq^s is set so that grain ro- 
tation is negligible on the time scale of solidification. 

Nucleation is incorporated into the simulations as fol- 
lows. Method I.: By including white noise into the gov- 
erning equations of amplitude that forces nucleation in 
the spatial and time windows used. Method II.: The sim- 
ulation area is divided into domains according to the lo- 
cal composition. The time-independent solution is found 
for these compositions. Critical fluctuations of statisti- 
cally correct numbers following Poisson distribution are 
placed into these areas in every time step. The added 
small-amplitude noise makes these critical fluctuations 
either grow or shrink. 

In nucleation-growth processes the transformed frac- 
tion often follows the Kolmogorov scaling, X(t) = 1 — 
exp{ — (t/to) p }, where the "Kolmogorov exponent" p is 
representative of the mechanism of the phase transfor- 
mation, and is evaluated from the slope of the plot 
ln[— ln(l — X)] vs. Int. In this work <j) < 0.5 is used 
to define the transformed fraction. 

First we apply the phase-field theory to predict nucle- 
ation rate in 3D. In the one- component Lennard- Jones 
system, the nucleation rate |1J and all the relevant phys- 
ical properties are known from molecular dynamics sim- 
ulations ||. The radial phase field profiles (Fig. 0a) 
indicate that the critical fluctuations are diffuse, and do 
not show bulk crystal properties for undercoolings larger 
than 14 K. The predicted interfacial free energy (Fig. |l]b) 
increases with temperature. While the phase-field predic- 
tions agree with results from computer simulations, those 
from the classical sharp interface theory differ from the 
experiments by eight to ten orders of magnitude. Similar 
results were obtained for the ice-water system (Fig. ||) 
with input data from ]l6[ . Without adjustable param- 
eters, a quantitative agreement has been achieved with 
computer simulations jl5| and experiment [|l7j, proving 
the power of the phase-field technique in attacking nu- 
cleation problems. 

In the case of binary alloys, such a rigorous test cannot 
be performed since the input information available for the 
crystal-liquid interface is far less reliable. In the nearly 
ideal Cu-Ni system, the critical undercoolings computed 
for the realistic range of nucleation rates (J = 10 -4 to 
1 drop _1 s _1 for electromagnetically levitated droplets of 
6 mm diameter) fall close to the experimental ones [Q 
(Fig. ||), indicating homogeneous nucleation. This con- 
tradicts the heterogeneous mechanism suggested earlier 



FIG. 2: Nucleation rate vs. temperature in the ice-water sys- 
tem. The experimental results (symbols) are from | ]17[ . The 
branches below 232 K indicate results obtained with physi- 
cally reasonable upper and lower estimates for the Gibbs free 
energy of undercooled water |2|. Notation as in Fig. M. 





FIG. 4: Soft-impingement in the phase-field theory, (a)- 
(d) 1000 x 1000 segments of 7000 x 7000 snapshots for the 
concentration field in cases (i)-(iv), respectively, taken at 
t/At = 4000, 1500, 1500, and 6000 (black and white corre- 
spond to the solidus and liquidus); (e) Kolmogorov exponent 
vs. crystalline fraction, n = X/Xm ax , where X max is the final 
value of the crystalline fraction; (f) experimental results for 
the crystallization of amorphous Fers.sSiir.sCuNbsBs M. 



FIG. 3: Nucleation temperature vs. composition that the 
phase field theory predicts for the nearly ideal Cu-Ni system. 
Upper and lower solid lines correspond to nucleation rates of 
10 -4 and 1 drop _1 s _1 for droplets of 6 mm diameter. The 
experimental data (squares) refer to electromagnetically lev- 
itated droplets pl| . The calculated liquidus and solidus lines 
(dashed) are also shown. 

]l8[ on the basis of Spaepen's value ans = 0.86 Q 
for the dimensionless interfacial free energy for the hard- 
sphere system. (For definition of a see |J.) Recent com- 
puter simulations |2(J] yielded considerably smaller values 
a HS = 0.51 and a^i = 0.58 (~ 0.6 we used), invalidating 
the earlier conclusion. These findings raise the possibility 
that homogeneous nucleation is more common in alloys 
than previously thought. 

We turn now to the problem of soft-impingement that 
we investigate in 2D using the properties of Cu-Ni alloys. 
Owing to the known difficulties of phase field simulations 
due to different time and length scales of the fields, we 
used an enhanced interface thickness (5 = 41.6 nm), a re- 
duced interfacial free energy (a = 0.1), and an increased 
diffusion coefficient (~ 100 x Dl). To ensure reasonable 
statistics and negligible influence of the periodic bound- 
ary conditions, the governing equations were solved nu- 
merically on a 7000 x 7000 grid under conditions (2lJ that 
ensure interfaces consisting of more than ten grid points 
needed for numerical accuracy. 

We begin by comparing three simulations where 
the particles were nucleated by Method II.: (i) large 
anisotropy so = 0.25 and small dimensionless nucleation 
rate J = 0.49 yielding 886 dendritic particles (the first 
large scale simulation of multiparticle dendritic solidifica- 
tion) [Fig. f|(a)]; (ii) sq = 0.25 with large nucleation rate 



J = 24.5 (10623 particles) [Fig. g(b)]; and (iii) isotropic 
growth with J = 24.5 (10528 particles) [Fig. 0(c)] . The 
respective Kolmogorov exponents differ from p = 2 that 
standard references Q assign for steady-state nucleation 
and diffusion controlled growth [Fig. 0(d)]. Our predic- 
tion for dendritic solidification, p « 3, obeys the relation- 
ship p = 1 + d (constant nucleation and growth rates in 
^-dimensions) confirmed experimentally (^3| , that follows 
from the steady-state traveling tip solution of the diffu- 
sion equation. In cases (ii) and (iii), p deviates from 2, as 
the formation of the diffusion layer is preceded by a tran- 
sient period in which phase field mobility controls growth. 
This period appears as an effective delay of the diffusion 
controlled process, yielding an initially enhanced p that 
decreases with increasing transformed fraction ||. This 
effect is pronounced at large nucleation rates for which 
the delay is comparable to the solidification time. In- 
deed, qualitatively similar behavior is observed at the 
extreme nucleation rates that occur during the forma- 
tion of nanocrystalline alloys [|J. 

The introduction of large amplitude noise into the gov- 
erning equations (Method I.) leads to comparable results. 
For s = 0.25 with C/>,o = 0.015, that yield - 830 den- 
dritic particles [case (iv)], more irregular shapes are pro- 
duced [Fig. 4(d)], and transient nucleation of "induction 
time" r/ At w 1400 is observed. The latter leads to fur- 
ther increased p, that reduces to ~ 3 when replacing t 
by t — r. Due to numerical stability problems appearing 
at large noise amplitudes, Method I. can be applied far 
from equilibrium, where nucleation occurs in reasonable 
simulation time and area. 

Summarizing, we demonstrated that the present phase 
field model is able to quantitatively describe crystal nu- 
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cleation in one- component 3D systems. The other predic- 
tions, including binary nucleation in 3D and transforma- 
tion kinetics in 2D, are also consistent with experiment. 
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